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Abstract 


The Aircraft Scheduling Problem (ASP) refers to allocating each aircraft to the optimal take-off and landing time and the 
appropriate runway. This problem is the allocation of aircraft to the desired runway so that the total damage due to delays 
or haste in landing or take-off of all aircraft is minimized. Runway allocation, landing and take-off sequences, and 
scheduling for each aircraft must be done in a predetermined time window. Time should also be considered as the time 
of separation between landings and take-offs due to the wake vortex phenomenon. In general, the purpose of such 
problems is to make maximum use of the runway. Therefore, in this study, a mathematical model of robust landing and 
take-off scheduling at an airport is provided, assuming no access to the airport runway at certain times. Moreover, delays 
and haste in landing and take-off on the runway, limited access to aircraft, runway repair time, and the possibility of 
runway disturbances are investigated. Robust optimization is used to deal with uncertainty at take-off and landing times. 
Finally, Genetic and Imperialistic Competitive Algorithm (ICA) are used to evaluate and analyze the problem because it 
is NP-HARD problem. The results indicate the ability of the proposed algorithms to find high-quality solutions in a 
short computation time for problems up to 7 runways and 60 aircraft. 


Keywords: Aircraft scheduling problem, Robust optimization, Genetic algorithm, Imperialistic competitive algorithm. 


1 | Introduction 
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Nowadays the high volume of transportation is done by aircraft (plane), which is very expensive 
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Air transportation plays a significant role in the movement of passengers and freight, and this role is 
becoming more important day by day so that air transportation accounts for a significant share of the 
transportation problem. There is a growing demand for air travel so itis not unexpected that the demand 
for air travel will multiply in the coming years. This has created many problems in recent years such as 
delays, traffic, and air accidents. 


The scheduling problem is one of the most important problems in the world today that has a significant 
impact on increasing the efficiency of production and transportation systems Zeng et al. [1]. Scheduling 
refers to the allocation of limited resources to activities that require that resource and is a kind of 
decision-making activity that is done to optimize one or more goals Ng et al. [2]. As a result, aircraft 
scheduling is one of the most important problems and factors affecting the transportation, industry and 
customer satisfaction. This problem is to determine the optimal sequence of landing or take-off and 
assigning them to different and same runways, taking into account the limitations of parking spaces so 
that the total cost of flight delays is minimized Lambelho et al. [3]. In other words, the Aircraft 
Scheduling Problem (ASP) deals with the scheduling of aircraft on parallel runways. Each aircraft has 
its weight penalty and can operate on a runway when ready Zhang et al. [4]. The actual operation time 
cannot be before the ready time and after the deadline. This time should be as close as possible to the 
target time of the aircraft Yang et al. [5]. It is known as the time window from the ready time to the 
deadline. It is impossible to deviate from the deadline. The aircraft in question is not assigned to the 
runway and the scheduling is known as infeasible if the aircraft operation time exceeds the deadline 
Shakibayifar et al. [6]. The target time can be predetermined, but the operation time depends entirely on 
the landing or take-off conditions. Thus, an aircraft may not be able to land or take-off at or near its 
target time Çiftçi and Ozkir [7]. Priority is given to landing rather than take-off, and priority is also given 
to large aircraft Zhang et al. [4]. Besides, the required separation time between the aircraft in question 
and other aircraft must be considered. Separation time is considered to prevent the risk of landing or 
take-off disturbances Sama et al. [8]. When landing or taking off, each aircraft produces turbulence that 
causes scheduling disturbances. This turbulence varies depending on the type of operation and the size 
of the aircraft Sama et al. [8]. Hence, the problems of scheduling in this area fall into two general 
categories: 1) airline problems including fleet, passenger, crew and flight crew allocation, and aircraft 
maintenance. 2) Problems related to take-off and landing restrictions on one runway to take-off and 
landing on different runways on the one hand and restrictions such as aircraft size, number of 
passengers, emergency landing, maintenance and failure of runways, weather conditions, etc. which take 
the problem out of its initial state and disrupt the initial schedule. Hence, a rescheduling based on the 
new conditions is required. Some studies focus on different heuristic algorithms and methods for 
expressing flight schedules due to the complexity of aircraft landing and take-off. Esmaeilidouki et al. 
[9] solved a problem of routing and transporting hazardous materials. They used a Genetic Algorithm 
(GA) to pay attention to the time of transportation and routing of materials. 


Simultaneously with other parts of the world, flight scheduling in Iran seems necessary due to the 
dynamics of distance scheduling due to the increase in flights and, consequently, the expansion of flight- 
related problems on the one hand and the impossibility of developing airport runways in the short term 
and the high cost of these projects on the other hand. Moreover, constraints such as accessible runway, 
adverse weather conditions, emergency landing, fuel shortages, security accidents, the unpredictability 
of flight during landing in general and for Iranian airports in particular, waste of time for cargo delivery, 
delays and crowds on departure, and the proximity of Iranian airports to major cities that cause problems 
for residents neat airports such as air pollution, noise pollution, etc. and cause the impossibility of flying 
at night or incurring additional costs for distance travel, etc. make it necessary to pay attention to flight 
schedules under conditions of uncertainty. 


Finally, in Iranian airports, including Tehran airport, problems such as a wide flight range can be 
mentioned, which has caused the average flight time to be unreasonably extended. On the other hand, 
problems such as aircraft idle times can be mentioned. The airport would have achieved better quality 
for both flight and customer service if there was a proper flight schedule. Thus, the main problem of 
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this study is to develop a model for scheduling the time interval between flights so that more benefits and 
lower costs ate achieved. In this study, an attempt has been made to provide a regular schedule taking into 
account the dynamic conditions at Tehran airport using scientific methods. In this way, improvements can 
be expected in various areas such as time interval scheduling, robust scheduling, attention to various 


components, time and cost savings, and so on. 


2 | Literature Review 


Airlines play an important role in transportation today, and air transportation has a major share of 
transportation Salehipour et al. [10]. As demand for air transportation grows rapidly, global air traffic is 
expected to double over the next 15 years. It is even expected to triple in some places. Many large airports 
are at high risk of air traffic delays or will be soon Zhou et al. [11]. The first idea to solve this problem is 
to increase the capacity of airports, however, this has also become a major problem Arkind [12]. In the 
current airport system, runways are considered the primary bottleneck at the airport Idris et al. [13]. Because 
the landing and take-off of all aircraft depend on the condition of the runways. So, it can be argued that a 
slight increase in runway efficiency significantly increases the efficiency of the airport. Many researchers 
have focused on the runway system, and one of the most important operational problems is the ASP. This 
problem tries to make the most of the runway by scheduling the aircraft landing and take-off as much as 
possible in a predetermined interval. 


With the passage of the Air Transport Liberalization Act in 1978, the way was paved for fundamental 
changes in the aviation industry Sinha [14] so that airlines could determine the flight route and ticket price 
of each route themselves. This created a competitive environment between airlines Sama et al. [15]. Today, 
according to forecasts, demand for aviation is expected to grow at an annual rate of 3% to 5%, despite the 
recent recession. Therefore, increased traffic causes high obstruction in the terminal areas and delays in 
the arrival of aircraft and creates long queues in cargo departure areas Fernandes et al. [16]. 


In this way, total air traffic is steadily increasing, and the number of commercial aircraft used will double 
over the next two decades Lieder and Stolletz [17]. Matching more flights is a significant challenge given 
the current level of transportation at busy airports. Runway capacity is often a limiting factor that creates 
plans to offer additional flights at the airport. Also, the allocation of the fleet has its costs and rules, which 
are mentioned below. As such, airlines are forced to use powerful tools to make informed decisions, reduce 
costs, and increase their share of existing demand to survive and establish themselves in a competitive 
market. One of these tools is the use of systematic scheduling approaches based on mathematical logic in 
different sectors of airlines. One of the main factors affecting the use of runways is the implementation of 
the minimum separation between aircraft landing, which is rooted in safety considerations Lee et al. [18]. 
Wave vortices are circulating air masses that are produced as a result of the flight. These vortices can be 
dangerous for one aircraft below another. The vortices generated by larger aircraft are stronger but they 
will not affect the safe operation of the aircraft if the aircraft is at a reasonable distance from another. 
Besides, if the lower aircraft is lighter than the upper, these vortices will have a greater impact on the lighter 
aircraft. Therefore, the minimum required distance between aircraft depends on the weight of the upper 
and lower aircraft. As a result, effective scheduling will help prevent lighter aircraft from landing 
immediately after a heavier aircraft lands or takes off at Zhou et al. [11]. 


Flight scheduling is a major part of an airline's operations because a variety of factors must be considered 
in developing a schedule, and flight operations must be highly safe and reliable. Developing a flight 
schedule or sub-schedule usually consists of a reciprocating process. This process is very time-consuming 
to reach an operational flight schedule and is usually done over six months Fernandes et al. [16]. In addition 
to the safety issues that are the responsibility of air traffic control, other stakeholders are very interested in 
how the landing is scheduled. The downside is the concern for airlines and airports. Airport operations, 
such as missions to patrol and carry passengers’ luggage, require careful scheduling, and delays in a landing 
may have detrimental effects on similar operations for subsequent aircraft. Airlines prefer schedules that 
minimize fuel costs Zhang et al. [4]. In studies by Allahverdi et al. [19] and Ball et al. [20], different problem- 


solving approaches and objective functions were reviewed in the literature related to flight schedules. In 
patticular, the main objective functions include minimizing delays and costs. Costs ate calculated based 
on deviations from the nominal schedule in terms of delays. Regarding the arrival schedule, Zhan et al. 
[21] predicted and assigned the sum of the differences between the landing times of each aircraft. Sama 
et al. [15] limited aircraft delays and deviations from the timetable by considering aviation priorities 
(aviation priorities refer to the classification of aircraft for a flight on the runway. For example, 
international aircraft are given priority due to the importance of on-time flights, and domestic flights are 
next). Khaksar and Sheikholeslami [22] presented a method for airline delay prediction via machine 
learning algorithms based on data mining, random forest Bayesian classification, K-means clustering, 
and hybrid approach and calculated delay occurrence and magnitude in both Iran and USA networks. 
In a study by Beasley et al. [23], the problem of aircraft displacement was investigated by considering 
transportation costs using the distance-to-cost adjustment approach. In particular, additional penalties 
should be considered if the aircraft is delayed concerning the initial solution. Bennel et al. [24] allocated 
landing times for aircraft based on the evaluation of a specific route by minimizing flight time and using 
routing patterns by maximizing the minimum time between two landings. Hu and Di Paolo [25] 
considered two objective functions: minimizing the delay of all aircraft and maximizing the length of all 
arrival queues. The first objective function emphasizes the operating cost of the aircraft, but the second 
objective function emphasizes the objective functions, which focus on the efficiency of using the airport 
capacity. Bennell et al. [26] examined the dynamic landing schedule at a single-runway airport. They 
considered the time window limit for each landing time and the minimum separation distance between 
successive landings where the separation time depends on the weight classes of the two aircraft at 
landing. The multi-objective formula was proposed. 


According to runway capacity, delays and fuel costs due to aircraft maneuvering, and additional flight 
time to reach the landing schedule. Kenan et al. [27] investigated the landing and take-off schedules of 
aircraft in terms of runway interdependencies and their heterogeneity and proposed an optimization 
method for the flight scheduling problem with general runway configurations (it is the runway that is 
used by 80% of flights, and there are some dedicated runways for private or military aircraft or giant 
aircraft). Cheng et al. [28] introduced the concept of flight operation risk assessment system for an airline 
and discussed the correlation between a risk factor and its sub components with fuzzy inference system. 
They also developed algorithms to identify the critical risk factors based on sensitivity of the risk factor 
and heuristic search. 


Kim et al. [29] developed a mathematical model to provide a basis for planning airport facilities and 
increasing existing services in Incheon International Airport. It was found that the proposed model has 
good prediction capability for traffic volume on the Incheon International Airport Expressway on an 
hourly basis. 


Tavakkoli Moghaddam et al. [30] considered the issue of aircraft meeting scheduling as a single band 
and solved it using a fuzzy planning approach. In their work, they proposed methods for the problem 
of aircraft landing with the shortest waiting time in the desired time window in critical situations, such 
as the closest landing time to the target time for each aircraft or the minimum landing time of aircraft. 


Rashidi Komijan et al. [31] present a mathematical model for an integrated airline fleet assignment and 
crew scheduling. They use Vibration Damping Optimization (VDO) algorithm to an appropriate 
solution to their problem during a reasonable time period. An experimental design based on the Taguchi 
method was taken into account too. In the discussion of solution methods. Hassanpour et al. [32] 
present a robust bi-level programming model for designing a closed-loop supply chain and use a robust 
bi-level and GA for solve the problem. 


The landing and take-off of each aircraft must be assigned to a runway or time. However, the 
requirements for separating two aircraft are met depending on the aircraft tail and minimizing the cost 
of delays. Some runways can only be used for landing, take-off, or certain types of aircraft. At the same 
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time, additional separation restrictions should be considered in scheduling airport runways. The dynamic 
scheduling approach presented in this study solves realistic examples of the problem for optimization in 
short computation time. Moreover, this study proposes a scheduling horizon for large examples that give 
near-optimal results. Therefore, in a comprehensive flight schedule, the modeling process is complex and 
difficult due to the large size of the problem, the types of constraints, the types of objective functions to 
be optimized, the parameters, and the types of decision variables. This is why researchers divide the 
problem into several smaller, independent problems and then solve each one separately. 


The main problems proposed in the flight schedule are as follows: 


I. Aircraft landing and take-off schedule design. 
Il. Fleet allocation to flight. 
III. Determining the flight path and aircraft maintenance schedule (allocation of aircraft to flight). 
IV. Flight crew scheduling. 


However, each of the more detailed problems needs to be optimized in each of the sub-problems. The 
flight scheduling problem, which covers all sub-problems, will not necessarily be optimal if the optimal 
solution for each sub-problem is discrete. In addition to the problem, in such a case, there is still a 
reciprocating process between each of the sub-problems, especially with the first stage, the landing and 
take-off schedule, so that the flight schedule is in a relatively optimal state, which will waste a lot of time. 
Therefore, adopting a different approach that considers all sub-problems as discrete can provide a more 
optimal solution to the flight scheduling problem and significantly reduce the time interval of the 
scheduling process. So, it seems necessary to adopt such an approach. 


In this study, a robust mathematical model is designed and solved for the problem of optimal scheduling 
of passenger aircraft flight time interval at Tehran airport. 


3 | Mathematical Modeling 


The ASP is to determine the optimal sequence of landing or take-off and assigning them to different and 
same runways, taking into account the limitations of parking spaces so that the total cost of flight delays is 
minimized. In other words, the ASP deals with the scheduling of aircraft on parallel runways. Hence, the 


modeling assumptions are as follows: 


I. The departure time of each aircraft depends on their sequence on the runway. 
Il. Each runway can direct a maximum of one aircraft at any moment. 
HI. Each aircraft can take-off or land on a maximum of one runway at any moment. 
IV. Runways ate not continuously available. 
V. All flight times are as robust uncertainty. 


Indices 

— i runway 

— j: aircraft (j) 

— I: aircraft (1) 

— t: time interval 
Input parameters 


n: the number of flights. 


m': the number of ready runways in each period t. 


dj: The flight time of aircraft j. 

Pi: The flight duration of aircraft | on the runway i. 

t; Access time for aircraft j to start the take-off process. 

a: Cost of haste in the landing/take-off of aircraft j. 

Bi: Cost of delays in landing/take-off of aircraft j. 

Mit: If it were possible for aircraft j to fly or land on runway i in period t, 1, and, otherwise, 0. 
Si: The ready time for aircraft | when ready to fly/land on one of the runways after aircraft j in step t. 
chi: The ready time for aircraft | when the first aircraft to take off or land on runway i in stage t. 
avf: The time to access runway i in period t. 

pd; : The probability of runway disturbances after landing/take-off of aircraft j in period t. 

Rt: The time required to repair the runway in each period t. 

M’: A large positive integer number. 

Decision variables 

Xj: If aircraft j flies/lands on runway i in period t, 1, and, otherwise, 0. 

Yir If aircraft j flies/lands on runway i in period t before aircraft 1, 1, and, otherwise, 0. 

Ci: The time to complete the flight of aircraft j in period t. 

Cj: The time to complete the flight of aircraft | in period t. 

Ej: The haste time of aircraft j for landing/take-off. 

Tj: The delay time of aircraft j for landing/take-off. 


The model variables including E, C, and T have non-negative integers, and variables X and Y have 


values of zero and one. 


3.1 | Decision Variables 


The objective function of the model is equal to the sum of the acceleration and delay times of the flights. 
In this formula, the values of the decision variables Ej and Tj for j-type aircraft are obtained from Eqs. 
(1) and (2), respectively. 
— t 
E; = max(0,d;- CĦ), (1) 
= t 
Tj = max (0,C} - dij), 2 
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se i ¥5 < Xj forallij,t j#1, 4) 
dito Yà = Xi foralltil 1#j, 5) 
Sh, YH =1 ~ forallt,i, 6) 
Xi < < < Mj for all tij, 7) 
CÌ- Ct> Si+ Pa» Yfi+ pdt +R'+(¥{- 1)*M forall t,,),1, 8) 
CÍ >0 forallt,j, 9) 
Ci- Crt > DE Dfa Yin Si + Dei chi* Yi + DE Pas V+ LPa 

Yj + De = 2p nent forallt,l j#1, 
Ce pn nav! ta Vint De ni cht e Y$ + Se nPa+Yi, forall tl, 11) 
C? > rj for allj, 12) 
Tj 2 Cj -dj forallt,j, 13) 
Tj 20 for all j, 14) 
Ej > dj- Cj for all t,j, 15) 
E; 20 for allj. 16) 


3.2 | The Limitations of the Proposed Model 


Constraint (3) indicates that each aircraft can land/take off on one and only one runway in each period t. 
Constraints (4) and (5) ensure that each aircraft has a landing/take-off operation immediately before and 
after only one other aircraft on the runway in each period t. Constraint (6) specifies the first flight on runway 
i in each period t. Constraint (7) introduces limited access to the runway. As stated in the model input 
parameters section, if aircraft j can fly/land on runway i in period t, the value of parameter Mi will be 1. 

Otherwise, it will be 0. The possibility of flight/landing of aircraft j on runway i in period t is e 
according to the movement set of aircraft j, i.e. Mi, which is a subset of runways M‘ and includes all runways 
on which aircraft j can land/take off. This constraint, thus, constrains the model to consider Mij which is 
one of the input parameters of the model, for assigning runway i to aircraft j and, consequently, assigning 
the value of 1 to the decision variable Xj. This assignment is made if the value of Mi is 1 such as X4. 
Constraint (8) ensures that the landing/take-off time of the aircraft is the result of the time it takes for the 
aircraft to be ready for flight, the duration of flight on the runway, and the likelihood of flight disturbances. 
Constraint (9) ensures that the flight time starts from zero. Constraint (10) ensures that the flight schedules 
do not interfere with each other. Constraint (11) ensures that the first aircraft does not have a problem 
flying. Constraint (12) specify the completion time of each flight. Constraints (13) and (14) specify the 
constraints on delay times for flight j. Constraints (15) and (16) specify the constraints on haste times for 
flight j. 


3.3 | Bertsimas and Sim Robust Uncertainty Approach 


The robust optimization specifies a suitable uncertainty set for imprecise input data and gives a solution 
that ensures feasibility in all amounts of uncertain parameters within the uncertainty set Ben-Tal and 
Nemirovski [33]. 


According to the evaluation, uncertainty in the flight time is added to the model using robust scheduling 
and the Bertsimas and Sim approach. With this change, Constraint (8) is modified as the Bertsimas model. 
Therefore, the proposed model is linear. Studies show that flight time is one of the important parameters 
whose values may exceed the nominal values. As a result, considering this parameter in uncertainty 
conditions can bring the proposed model closer to the reality of the problem. Robust scheduling and the 
Bertsimas and Sim approach are used to consider uncertainty in the flight time. Robust optimization seeks 


optimal or near-optimal solutions that are likely to be feasible. The Bertsimas and Sim approach is one 
of the four main approaches for considering uncertainty in robust scheduling. This section briefly 
mentions this apptoach. For this purpose, the following linear programming model is considered. 
Min jj cj Xj, 
s.t. (17) 
Ax <b. 


In this model, it is assumed that only the right-hand coefficients in the constraints, matrix A, have 
uncertain values, and the values of this matrix, i.e. g _, fluctuate in the range / dy -â ,ã +â], where 4. 
i 5 8g Ay ta, j 
and 4 are the nominal values and the maximum deviation of parameter a, , respectively. The proposed 
ij 

Bertsimas and Sim robust model is as follows: 

min ) CjX;, 

j 
s.t. 


1 
X; ajx + zi Ty + Mig, by <b; for alli, (18) 
j 45 i jJi Hij 


Zit Hij 


Zi; Hij 


aj Xij for all i, j, 
0 for alli,j, 


iv IV 


where z, and //; are dual auxiliary variables, and parameter r , called the uncertainty budget, indicates 


the level of conservatism that is chosen according to the importance of the constraint as well as the risk- 
taking of the decision-maker. Hence, Constraints (8), (10) and (11) become Constraints (19) to (22) 
according to the Bertsimas and Sim robust model. 

Pa: The tolerance of flight time of aircraft 1 on runway i. 

r; The uncertainty budget for flight time. 

Variables 

ppi and qin: The robust model variables. 


Ci- Cf> Si+ Pu» Yh + Tuppa + que + pdf «R'+(Yf - 1)« 


19) 
M for all t,i,j,1, 
= t t t 
Ci- Ch > DR DE Yit Si + L chi « Yio + De Dia Pa * Yi + 20) 
t t è 
Em Pi t Yi + Ppp + qii + Xizi Lia pd; «Rts Yj for allt,l j #1, 
t it t 
Cf > DEravi* Yi + DE chi + Yi + DE Pa» Yio + DE Pa* Yio +l pp + 21) 
gut for allt,1, 
Ppa + dit > Pu* Yi for alli,Lt. 22) 


4 | Funding 


In this section, the proposed mathematical model is validated. The collected data indicate that the 
movement position of the aircraft does not follow the same distribution functions. The status of the 
mathematical model parameters is as follows based on the timing performed. 
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Distribution Summary 


Distribution: Triangular 
expression: TRIA(10, 12.5, 18) 
Square Error: 0.000380 


Fig. 1. Aircraft ready time on the runway. 


The ready time parameter of aircraft } when selected as the first aircraft. 


Distribution Summary 


Distribution: Triangular 
Expression: TRIA(75, 90.2, 117) 
Square Error: 0.000137 


Fig. 2. Aircraft's ready time when selected as the first aircraft. 


According to the evaluation, the aircraft's ready time for landing and flight follows the triangular function 
with a minimum value of 10 minutes and a maximum value of 18 minutes, which is coded based on the 
obtained values of the mathematical model. 


The ready time parameter of aircraft j when selected as the first aircraft. 


Table 1. Adjustment of the mathematical model parameters. 


The Mathematical Model Parameters Intended Values 
The flight time of aircraft j U~ [10,30] 

The flight time of aircraft j on the runway U= 3] 

The time to access aircraft j to start the take-off process U~ [1,6] 

The cost of hastening the landing/take-off U~ [20,50] $ 

The cost of Delays in landing/ take-off U~ [80,100] $ 


The possibility of runway disturbances after landing/takeoff 0.35 
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According to the adjustment of the mathematical model parameters, the flight schedule of one of the 
airports is considered as follows (it is worth mentioning that small, medium, and large runways are selected 
for the schedule. Therefore, 20 aircraft are investigated in 3 runways over a 24-hour schedule. The flight 
schedule is as follows: 


Small runway 


Medium 
runway 


Large runway 


14 15 16 17 18 19 20 21 22 23 24 
Medium 
runway 


Fig. 3. The flight schedule for 24 hours. 


The evaluation of the designed problem and flight time of the aircraft show that 50% of the aircraft fly 
on small runways, 30% on medium runways, and 20% on large runways. The flight schedule for 24 
hours is shown in Fzg 3. Since the problem under study is NP-HARD, in this section, the mathematical 
model is evaluated using genetic meta-heuristic and imperialist competitive algorithms. 


4.2 | Adjusting the Meta-Heuristic Algorithm Parameters 


One of the most important steps in designing meta-heuristic algorithms to achieve optimal solutions is 
algorithm calibration. Different values of control parameters of these algorithms may affect computation 
indices including solution quality and computation time, so a series of calibration tests are often 
performed to find the optimal combination of different values of algorithm control parameters. The 
algorithms proposed in this study are genetic and imperialist competitive algorithms. GA control 
parameters include initial population size (Popga), crossover rate (P,), mutation rate (P,,), and the 


maximum number of generations (G The imperialist competitive algorithm control parameters are 


ax) 
the number of neighborhoods produced (Pop;c 4), the number of iterations (Decade), and the percentage 
effect of the total power of each imperial colony on its power (PICA). Each of these parameters affects 
the computation indices in a certain range of their values and has a negligible effect outside this range. 
The computation indices used in the tests of this section and future sections are the mean values of the 
objective function for the optimal solutions and their mean computation times per ten iterations of the 
algorithm. A set of tests based on the Taguchi method is designed to investigate the interaction of 
control parameters of the proposed algorithms and achieve their optimal combination. Each of the 
parameters for each of the proposed algorithms in the previous section is tested at three levels. Solution 
levels are the mean solutions obtained and the mean computation times. To simultaneously consider the 
quality of solutions and computation times, the values of the two are normalized and added together. 
First the parameters of the GA and then the parameters of the Imperialistic Competitive Algorithm 
(ICA) algorithm were examined. Tab/e 2 shows the desired factors and their levels for the GA. The tests 
required to examine the various combinations of factors and the corresponding solutions can be seen 
in Table 3. The generated data are analyzed by MINITAB 14 software, and the results are provided in 
the following tables. 


Table 2. Factors and their levels. 


Factor Levels 
Popga 200, 300, 400 
Po 0.8, 0.7, 0.6 
Pa 0.15, 0.12, 0.10 
Grax 200, 300, 400 
Pii 0.10, 0.15, 0.2 


To simultaneously examine the effect of factors on the quality of solutions and computation times, their 
values are normalized, added together, and become a solution variable. Finally, the inverse of this value 
is calculated and considered as the solution variable. The larger the value, the better. This can be seen 
in Table 4. 
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Table 3. Combinations of factors and corresponding solution levels in multifactorial tests. 


IJRIE Popca Pe Pm Pma Gmax Solution CPU(s) 
200 0.6 0.10 0.10 200 164.00 6.56 
200 0.6 0.10 0.10 300 42.00 11.80 
200 0.6 0.10 0.10 400 46.70 13.73 
31 200 0.7 0.12 0.15 200 40.00 8.50 
200 0.7 0.12 0.15 300 13.67 12.60 
200 0.7 0.12 0.15 400 32.33 13.75 
200 0.8 0.15 0.20 200 66.33 12.45 
200 0.8 0.15 0.20 300 26.33 16.58 
200 0.8 0.15 0.20 400 27.00 20.70 
300 0.6 012 0.20 200 20.67 14.17 

300 0.6 012 0.20 300 14.33 23.07 

300 0.6 012 0.20 400 48.33 12.80 

300 0.7 015 0.10 200 20.33 16.21 
300 0.7 015 0.10 300 21.67 22.89 

300 0.7 0.15 0.10 400 44.67 13.21 

300 0.8 0.10 0.15 200 18.67 22.58 

300 0.8 010 0.15 300 22.33 23.54 

300 0.8 010 015 400 49.10 12.51 

400 0.6 015 0.15 200 15.00 28.12 
400 0.6 015 0.15 300 45.20 14.30 
400 0.6 015 0.15 400 26.57 20.01 
400 0.7 0.10 0.20 200 14.87 27.95 
400 0.7 0.10 0.20 300 14.65 27.80 
400 0.7 0.10 0.20 400 13.00 29.14 
400 0.8 0.12 0.10 200 8.60 34.00 
400 0.8 0.12 0.10 300 9.20 31.23 
400 0.8 0.12 0.10 400 7.60 31.46 


Table 4. Combinations of factors and normalized solution levels in multifactorial tests. 


Popea Pe Pm Pmu Gmax Normalized 1/Normalized 
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200 0.6 0.10 0.10 200 0.200401 4.99 

200 0.6 0.10 0.10 300 0.070722 14.14 
200 0.6 0.10 0.10 400 0.079804 12.53 
200 0.7 0.12 0.15 200 0.062105 16.10 
200 0.7 012 0.15 300 0.039810 25.12 
200 0.7 0.12 0.15 400 0.063385 15.78 
200 0.8 0.15 0.20 200 0.099832 10.02 
200 0.8 0.15 0.20 300 0.061938 16.15 
200 0.8 0.15 0.20 400 0.070604 14.16 
300 0.6 0.12 0.20 200 0.050836 19.67 
300 0.6 0.12 0.20 300 0.060636 16.49 
300 0.6 0.12 0.20 400 0.079888 12.52 
300 0.7 0.15 0.10 200 0.054357 18.40 
300 0.7 0.15 0.10 300 0.068697 14.56 
300 0.7 0.15 0.10 400 0.076483 13.07 
300 0.8 0.10 0.15 200 0.064667 15.46 
300 0.8 0.10 0.15 300 0.070699 14.14 
300 0.8 0.10 0.15 400 0.080214 12.47 
400 0.6 0.15 0.15 200 0.071084 14.07 
400 0.6 0.15 0.15 300 0.079179 12.63 
400 0.6 0.15 0.15 400 0.068788 14.54 
400 0.7 0.10 0.20 200 0.070609 14.16 
400 0.7 0.10 0.20 300 0.070070 14.27 
400 0.7 0.10 0.20 400 0.070749 14.13 
400 0.8 0.12 0.10 200 0.075026 13.33 
400 0.8 0.12 0.10 300 0.070403 14.20 


400 0.8 0.12 0.10 400 0.069012 14.49 


Table 5. Estimated correlation coefficients of the model for SN ratios. 


Term Coef SE Coef T P 

Constant 22.9674 0.4156 55.270 0.000 
Popga200 -0.4973 0.5877 -0.846 0.410 
Popga300 0.5624 0.5877 0.957 0.353 


Pe 0.6 -0.8177 0.5877 -1.391 0.043 
Pe 0.7 1.0433 0.5877 1.775 0.025 
Pam 0.10 -1.1060 0.5877 -1.882 0.038 
P,0.12 1.1476 0.5877 1.953 0.079 
Pry 0.12 -0.9031 0.5877 -1.537 0.144 
Pamu 0.15 0.7053 0.5877 1.200 0.248 
Gmax200 -0.5723 0.5877 -0.974 0.345 
Gmax300 0.8050 0.5877 1.370 0.009 


Factor correlation coefficients for SN ratios are presented in Tabe 5. Coefficients with larger absolute 
values are more important than other factors. These coefficients are used to rank the factors in the 
solution table. As can be seen in the table, the factors P, = 0.6, P,, = 0.10, P pu = 0.12, and Gmax = 
300 have a significant effect on the solutions at the 95% confidence level. 


Table 6. Analysis of variance for SN ratios. 
Source DF SeqSS AdjSS AdjMS F P 


FOP size 2 5110 5.110 2.555 0.55 0.589 
P. 2 16.272 16272 8.136 1.75 0.016 
Pin 2 22877 22.877 11.439 2.45 0.018 
P mu 2 12.170 12.170 6.085 1.31 0.299 
Ce. 2 9267 9.267 4.633 0.99 0.029 


Residual Error 14 74.599 74.599 4.662 
Total 24 140.296 


Similar results are presented in Tab/e 7, which shows the factor correlation coefficients for the mean 
solutions. Moreover, the analysis of variance for the SN coefficient and the mean solutions are 
performed for the tested factors, and the results are given in Tables 6 and 8. According to the analysis of 
variance, the factors P,,, P, and G,,,;with a p-value of less than 0.05 have a significant effect on the 
solutions as expected. 


Table 7. Estimated correlation coefficients of the model for mean solutions. 


Term Coef SE Coef T P 
Constant 14.5035 0.3741 38.765 0.000 
Popga200 -0.1711 0.8344 -0.205 0.840 


Popga300 0.6944 0.8344 0.832 0.417 
P. 0.6 -0.9944 0.8344 -1.192 0.251 
P0.7 1.6733 0.8344 2.006 0.042 
Pa 0.10 -1.5822 0.8344 -1.896 0.016 
P0.12 1.9078 0.8344 2.287 0.036 
Pim 01 -1.2022 0.8344 -1.441 0.169 
Pim015 1.0867 0.8344 1.302 0.211 
Gmax 200 -0.4811 0.8344 -0.577 0.042 
Gmax300 1.2411 0.8344 1.487 0.000 


Table 8. Analysis of variance for mean solutions. 


Source DF SeqSS AdjSS AdjMS F P 
Popga 2 7.069 7.069 3.534 0.38 0.692 
: 2 38.249 38. 19.124 2.03 0,022 
Pin 2 56.241 56.241 28121 2.99 0.006 
Penis 2 23.756 23.756 11.878 1.26 0.309 
ome 2 21145 21145 10572 112 0.001 


Residual Error 14 150.370 150.370 9.398 
Total 24 296.829 
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Table 9. SN ratio solution. 


Level POPsize P: Pin Piru Gmax 
1 22.47 22.15 21.86 22.06 22.40 
2 23:53 24.01 24.12 23.67 23.77 
3 22.90 22.74 22.93 23.17 22.73 
Delta 1.06 1.86 2.25 1.61 1.38 
Rank 5 2 1 3 4 


Solution levels are evaluated according to mean solution indices and SN ratios, and factors are ranked to 


determine the priority or degree of importance of each factor. The mean solution index for each level of 


each factor can be seen in the solutions table. The ranking of the factors according to the solution analysis 


concerning SN coefficients and means can be seen in Tabes 9 and 10. Accordingly, the P,,, factor has the 
highest rank in both tables, and the P, factor has the second rank. The ranking of other factors is the same 


for the mean solution indices and SN coefficients. 


Table 10. Mean solutions. 


Level PODgize Pe Pin Pinu Gmax 
1 14.33 13.51 12.92 13.30 14.02 
2 15.20 16.18 16.41 15.59 15.74 
3 13.98 13.82 14.18 14.62 13.74 
Delta 1.22 2.67 3.49 2.29 2.00 
Rank 5 2 1 3 4 


Factor interaction analysis is used to fine-tune the levels of factors. These effects can be seen in the 


following figures (Figs. 4 and 5). 


Mean of Means 


Main Effects Plot (data means) for Means 
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Fig. 4. Mean solutions. 


As shown in Fig. 4, the mean solution index is maximized by Pope, factor at level 300, by P, factor at level 
0.6, by P „factor at level 0.12, by Pnu factor at level 0.15, and by Gmax 
levels of the factors in Fg. 5 shows that the factors maximize the SN index at similar levels. 


Mean of SNratios 


factor at level 300. Examining the 


Main Effects Plot (data means) for SN ratios 
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Fig. 5. Mean SN ratio. 


Thus, the values of the control parameters of the proposed GA are adjusted according to the following 


table. 


Table 12 shows the factors and their levels for the imperialist competitive algorithm. The tests required 
to examine the various combinations of factors and the corresponding solutions are presented in Tab% 


13. 


To simultaneously examine the effect of factors on the quality of solutions and computation times, their 
values are normalized, added together, and become a solution variable. Finally, the inverse of this value 
is calculated and considered as the solution variable. The larger the value, the better. This can be seen 


in Table 14. 


Table 11. Optimal values of control parameters. 


Table 13. Combinations of factors and solution levels in multifactorial tests. 


Factor Levels 
Popga 300 

P: 0.7 

Pi 0.12 
Pr 0.15 
Giris 300 


Table 12. Factors and their levels. 


Factor Levels 
Popca 300, 400, 500 
decade 200, 300, 400 
Pica 0.08, 0.10, 0.12 


Popica Decade Pica Solution CPU (s) 
300 200 0.08 85.2 0.66 
300 200 0.08 40.6 0.82 
300 200 0.08 40.2 1:13 
300 300 0.10 36.0 0.65 
300 300 0.10 23.0 0.83 
300 300 0.10 27.8 1.14 
300 400 0.12 46.2 0.65 
300 400 0.12 31.4 0.80 
300 400 0.12 30.8 1.15 
400 200 0.10 19.0 1.29 
400 200 0.10 19.0 1.86 
400 200 0.10 25.6 0.98 
400 300 0.12 30.6 1:27 
400 300 0.12 37.8 1.87 
400 300 0.12 36.2 0.99 
400 400 0.08 19.0 1.28 
400 400 0.08 21.8 1.92 
400 400 0.08 23.0 0.98 
500 200 0.12 21.4 2.55 
500 200 0.12 24.6 1.25 
500 200 0.12 19.0 1.77 
500 300 0.08 22.6 2.59 
500 300 0.08 26.0 1.24 
500 300 0.08 23.4 1.66 
500 400 0.10 32.2 2.60 
500 400 0.10 28.2 1.30 
500 400 0.10 22.6 1.74 
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Table 14. Combinations of factors and normalized solution levels in multifactorial tests. 


I] RIE Popica Decade Pica Normalized 1/Normalized 
300 200 0.08 0.122624 8.16 

300 200 0.08 0.072106 13.87 
300 200 0.08 — 0.080000 12.50 

35 300 300 0.10 0.061851 16.17 
300 300 0.10 0.050734 19.71 
300 300 0.10 0.065022 15.38 
300 400 0.12 0.074394 13.44 
300 400 0.12 0.060252 16.60 
300 400 0.12 0.068981 14.50 
400 200 0.10 0.058258 17.17 
400 200 0.10 0.073676 13.57 
400 200 0.10 0.057989 17.24 
400 300 0.12 0.071981 13.89 
400 300 0.12 0.097065 10.30 
400 300 0.12 0.071294 14.03 
400 400 0.08 0.057987 17.25 
400 400 0.08 0.078742 12.70 
400 400 0.08 0.054791 18.25 
500 200 0.12 0.095291 10.49 
500 200 0.12 0.064062 15.61 
500 200 0.12 0.071241 14.04 
500 300 0.08 0.097848 10.22 
500 300 0.08 0.065513 15.26 
500 300 0.08 0.073676 13.57 
500 400 0.10 0.109924 9.10 
500 400 0.10 0.069841 14.32 
500 400 0.10 0.074857 13.36 


Table 15. Estimated correlation coefficients of the model for SN ratios. 


Term Coef SE Coef T P 

Constant 22.5770 0.7015 32.184 0.000 
Popyca300 0.2102 0.9921 0.212 0.852 
Popyca400 0.5968 0.9921 0.602 0.009 


decade200 -0.3813 0.9921 -0.384 0.738 
decade300 0.1738 0.9921 0.175 0.045 
Pica0.08 -0.5050 0.9921 -0.509 0.661 


Pica0.10 0.6425 0.9921 0.648 0.037 


Table 16. Analysis of variance for SN ratios. 


Source DF SeqSS AdjSS AdjMS F P 
Popa 2 13.1550 13.1550 11.5775 4.79 0.026 
Popica 2 10.6559 10.6559 10.3279 4.22 0.031 
Pica 2 22.0603 22.0603 19.323 8.49 0.04 


Factor correlation coefficients for SN ratios are presented in Tab 15. Coefficients with larger absolute 
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values are more important than other factors. These coefficients are used to rank the factors in the solution 
table. As can be seen in the table, the factors Popjc, = 400, decade = 300, and Pic, = 0.1 have a significant 
effect on the solutions at the 95% confidence level. 


Similar results are presented in Table 17, which shows the factor correlation coefficients for the mean 
solutions. Besides, the analysis of variance for the SN coefficient and the mean solutions are performed 
for the tested factors, and the results are given in Tables 16 and 18. According to the analysis of variance, 
the factors, Pica, decade, and Pop;c with a p-value of less than 0.05 have a significant effect on the solutions 
as expected. 


Term Coef SE Coef T P 

Constant 14.0996 0.2973 47.433 0.000 
Popyca300 0.3801 0.4204 0.904 0.381 
Popyca400 1.0132 0.4204 2.410 0.030 
decade200 -0.4721 0.4204 -1.123 0.280 
decade300 -1.9308 0.4204 -4.593 0.000 
Pica0.08 -0.5689 0.4204 -1.353 0.197 
Pica0.10 1.4109 0.4204 3.356 0.005 


Table 18. Analysis of variance for mean solutions. 


Table 17. Estimated correlation coefficients of the model for mean solutions. 


Source DF SeqSS AdjSS AdjMS F P 
Popica 2 20.817 20.817 10.408 4.36 0.034 
decade 2 57.503 57.503 28.751 12.05 0.001 
Pica 2 53.902 53.902 26.951 11.30 0.001 
Residual Error 14 33.400 33.400 2.386 

Total 20 165.622 


Table 19. SN ratio solution. 


Level Popjca Decade Pica 
1 22.79 22.20 22.07 
2 23.17 22.75 23.17 
3 21.77 22.78 21.77 
Delta 1.40 0.59 1.15 
Rank 1 3 2 


Table 20. Mean solutions. 


Level Popjca Decade Pica 
1 14.48 13.63 13.53 
2 14.93 14.28 15.11 
3 12.89 14.39 13.66 
Delta 2.05 0.76 1.58 
Rank 1 3 2 


Solution levels are evaluated according to mean solution indices and SN ratios, and factors are ranked 
to determine the priority or degree of importance of each factor. The mean solution index for each level 
of each factor can be seen in the solutions table. The ranking of the factors according to the solution 
analysis concerning SN coefficients and means can be seen in Tables 19 and 20. Accordingly, the Popic, 
factor has the highest rank in both tables, and the Pjc, factor has the second rank. The ranking of 
decade factors is the same for the mean solution indices and SN coefficients. Factor interaction analysis 


is used to fine-tune the levels of factors. These effects can be seen in the following figures (Figs. 6 and 


7). 


As shown in Fig. 6, the mean solution index is maximized by Pop;c,factor at level 400, by decade factor 
at level 200, and by Pc, factor at level 0.1. Examining the levels of the factors in Fg. 7 shows that the 
factors maximize the SN index at similar levels. 
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Main Effects Plot (data means) for Means 
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Fig. 6. Mean solutions. 
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Fig. 7. Mean SN ratio. 


Thus, the values of the control parameters of the proposed imperialist competitive algorithm are adjusted 
according to the following table. 


Table 21. Optimal values of control parameters. 


Factor Levels 
Popica 400 
decade 300 


4.3 | Evaluation of Algorithms 


A Relative Percent Deviation (RPD) is used as a measure for comparing GA and ICA algorithms Eq. (23). 
Small, medium and large problems are used to measure the performance of these algorithms. The test 
results, the best solution, the worst solution, the mean solutions, and the mean execution time are presented 
in Table 22. Moreover, Table 23 shows the RPD values and the mean execution time. The condition for 
stopping the algorithms is considered to be 100 seconds. Preliminary tests show that these algorithms 
usually reach the best solution before this time, and the first time they reach the best value of the objective 
function is recorded. 


SOlavg—SOlmin 


RPD = (23) 


solmin 


Table 22. Values obtained from different executions for GAMs and both proposed algorithms. 


AP Runway GAMS GA ICA 
Global Time Best Worst Mean Time Best Worst Mean Time 
Optimal 
3 1 0:00:06 1 1 1 0:00:16 1 1 1 0:00:14 
4 5 4 0:00:09 4 4 4 0:00:17 4 4 4 0:00:17 
7 9 0:00:11 9 9 9 0:00:20 9 9 9 0:00:19 
3 6 0:00:10 6 6 6 0:00:30 6 6 6 0:00:30 
6 5 12 0:05:42 12 12 12 0:00:37 11 12 11.33 0:00:33 
7 11 0:19:11. 11 11 11 0:00:43 11 11 11 0:00:43 
3 15 0:56:01 15 15 15 0:01:21 15 15 15 0:01:47 
10 5 - - 20 26 21.2 0:02:33 21 26 22.2 0:02:29 
7 - - 18 22 19.4 0:03:00 14 17 15.6 0:03:21 
3 - - 47 51 47.66 0:21:36 48 51 49.8 0:23:04 
20 5 - - 46 48 46.2 0:29:58 49 50 49.5 0:27:53 
T - - 19 25 22 0:28:46 19 24 23.33 0:28:49 
3 - - 35 40 36.2 0:32:16 -35 39 37.8 0:30:10 
40 5 - - 34 37 35.6 0:38:36 36 40 39.6 0:39:23 
7 - - 16 18 17.2 0:36:34 15 17 16.8 0:37:12 
3 - - 70 86 72.2 0:37:45 69 73 71.2 0:36:19 
60 5 - - 99 61 56.6 0:38:44 58 64 63.6 0:40:55 
7 z - 68 73 70 0:40:08 68 74 70.2 0:41:17 
Table 23. RPD values and average execution times calculated. 
GA ICA 
AP Runway RPD Average Computation Time (Sec) RPD Average Computation Time (Sec) 
3 0.00 0:00:16 0.00 0:00:14 
4 5 0:00 0:00:17 0:00 0:00:17 
7 0:00 0:00:20 0:00 0:00:19 
3 0:00 0:00:30 0:00 0:00:30 
6 5 0:00 0:00:37 0:00 0:00:33 
7 0:00 0:00:43 0:00 0:00:43 
3 0:00 0:01:21 0:00 0:01:47 
10 5 0:06 0:02:33 0:06 0:02:29 
7 0:08 0:03:00 0:11 0:03:21 
3 0:01 0:21:36 0:04 0:23:04 
20 5 0:004 0:29:58 0:01 0:27:53 
7 0:16 0:28:46 0:23 0:28:49 
3 0:03 0:32:16 0:03 0:30:10 
40 5 0:05 0:38:36 0:10 0:39:23 
7 0:08 0:36:34 0:12 0:37:12 
3 0:03 0:37:45 0:03 0:36:19 
60 5 0:07 0:38:44 0:10 0:40:55 
7 0:03 0:40:08 0:03 0:41:17 


Statistical results indicate that GA performs better than the ICA algorithm. The mean plot with LSD 
intervals for the two algorithms is shown in Fz. 8 for further analysis of the results. According to the 
tables above and the figure below, the values of the objective function and the execution time in both 
GA and ICA algorithms are close to each other. The GA algorithm performs better than the ICA 
algorithm in terms of both the value of the objective function and the execution time. 
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Interval Plot of GA; ICA 
95% CI for the Mean 


Fig. 8. Mean plot and LSD intervals (95% confidence level) for GA and ICA algorithms. 


RPD values for different numbers of aircraft and runways are shown in Figs. 9 and 70. 
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Fig. 9. RPD values for different number of aircraft. 


0/6 
0/5 
0/4 


0 T T T 1 
0 2 4 6 


Different number of runways 


lee) 


Fig. 10. RPD values for different number of runways. 


As shown in Figs. 9 and 70, as the number of aircraft and runways increases, the difference between the 
performances of the two proposed algorithms is greatly reduced. 


5 | Conclusion 


This is conducted to solve the aircraft take-off and landing problem in situations where take-off and landing 
times are uncertain and aircraft also have access restrictions. It also considers the possibility of runway 
disturbances and the time required to repair the runway. None of the previous studies have examined the 
assumption of the possibility of runway disturbances, the time required to repair the runway, the lack of 
access to aircraft at certain times, and the uncertainty of take-off and landing times at the same time. On 
the other hand, finding the sequence of landing and take-off on runways to minimize the total delay and 


haste in take-off and landing of aircraft is one of the most difficult problems of combined optimization 
even on a medium scale. In this study, a new integer scheduling model for the flexible airline scheduling 
problem is provided to minimize the total weight of aircraft delays and haste considering constraints 
such as aircraft's ready time, runway disruptions, limited aircraft access, and runway access time. The 
computational difficulties of this problem have made it impossible to solve it accurately on a large scale. 
According to studies, this is an NP-hard problem even for airport scheduling mode, which becomes 
almost impossible to solve through precise methods and optimization software as the number of aircraft 
increases. So, the use of heuristic and meta-heuristic algorithms can be helpful. GA and ICA algorithms 
are used to solve the problem. 


We first designed the corresponding mathematical programming model to solve the problem studied in 
this research work. Then, we presented the robust counterpart of the problem model using a robust 
planning approach. We solved the problem with certain and uncertain conditions and discussed several 
scenarios for parameters involved in the solution. Finally, we solved the problem using genetic and 
imperialist competitive algorithms concerning its complicated structure. For this purpose, we proposed 
the GA and imperialist competitive algorithm to solve this problem and employed different parameters 
for them. 


Based on obtained results, the mathematical programming model can be solved in a shorter time and 
with sufficient accuracy. It is because that the data from the studied sample is such that we are not 
dealing with a large-scale problem at first, and metaheuristic algorithms cannot compete with the exact 
methods such as the branch and bound algorithm in terms of execution time and solution quality. 
Therefore, it is observed that the exact method solution is more accurate than that of the genetic and 
imperialist competitive algorithm employed in this study. But the results show the ability of algorithms 
which are applied to solve this model. The presented algorithm can provide acceptable solutions for 
problems, which cannot be solved by exact methods and different solvers such as CPLEX in a 
reasonable time. 


In this study, the problem is investigated in terms of the number of aircraft and runways to evaluate the 
algorithms used. The results confirm the correct operation of the algorithms. According to the results, 
further studies are recommended to use accurate problem-solving methods such as the Benders 
decomposition algorithm or variable neighborhood descent meta-heuristic algorithm. It is also 
recommended that fuzzy theory and logic of gray numbers be used in the mathematical model and that 
the results be compared with this study. 
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